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Abstract 

A numerical scheme is described for accurately accommodating oblique, non-aligned, 
boundaries, on a three-dimensional cartesian grid. The scheme gives second-order ac- 
curacy in the solution for potential of Poisson's equation using compact difference 
stencils involving only nearest neighbors. Implementation for general "Robin" bound- 
ary conditions and for boundaries between media of different dielectric constant for 
arbitrary-shaped regions is described in detail. The scheme also provides for the inter- 
polation of field (potential gradient) which, despite first-order peak errors immediately 
adjacent to the boundaries, has overall second order accuracy, and thus provides with 
good accuracy what is required in particle-in-cell codes: the force. Numerical tests on 
the implementation confirm the scalings and the accuracy. 

1 Introduction 

For certain types of modeling problem it is advantageous to use cartesian finite-differences 
on a lattice that does not conform to the interfaces or boundaries of the space under con- 
sideration. Instead of the common adoption of an unstructured mesh that conforms to the 
geometry, the interfaces are considered to be arbitrary surface shapes, and the difference 
equations are modified adjacent to them to account for their boundary conditions. Such 
a choice is natural, for example, if the interfaces actually move in time, since then the 
remeshing cost might become excessive. 

One option for cartesian finite difference is to regard the interfaces as being approximated 
by "staircase" contours that follow cell boundaries. The problem with such a choice is that it 
constitutes a gross approximation for smooth, possibly curved, diagonal interfaces. In view 
of the rapidly increasing computational cost of reducing the cell size in multiple dimensions 
so as to reduce the errors in such a representation, a more attractive option appears to be 
to accommodate diagonal interfaces with a difference scheme and appropriate interpolation 
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that is higher order, notably quadratic, in the cell size, but depends as far as possible only 
on the local values of the potential. 

This paper describes a scheme for solving Poisson's equation with a wide range of general 
boundary conditions at interfaces, and subsequently interpolating the potential and field 
(potential-gradient) back to arbitrary points in the solution region. The scheme has been 
built into a 3-D Particle-in-Cell (PIC) code for solving self-consistent electrostatic plasma 
problems[l|. Examples to demonstrate the convergence have been studied. 

2 One-dimensional differences and interpolation 
2.1 Difference approximations 

First consider a one-dimensional mesh (not necessarily equally spaced) of nodes Xi on 
which a potential quantity is defined. For cells not affected by interfaces, we consider 
the gradient to be given at positions half-way between nodes, Xj+1/2 = (a^i + Xi+i)/2, as 
— ~ (pi)/dxi, where dxi = Xj+i — Xi and to be interpolated linearly between 

adjacent values. See Fig. [TJ So for cells that are not affected by interfaces, (f)'{x) = 
Wi+i/2{^ - Xi-1/2) + 4>'i-i/2{xi+i/2 - x)]/{xi+i/2 - Xi^i/2). The second derivative is then 
'Pi = (0i+i/2 ~ 0i-i/2)2/(c^a;j + dxi-i), which we note is uniform because this is a second- 
order interpolation of 0. Because 0" is uniform, it introduces no inconsistency to regard the 
0" registration as being at position Xi, the same place as 0i, even if the mesh is non-uniform 
so that Xi is not half way between Xi-1/2 and Xi+i/2- 




X 



Figure 1: Mesh notation in the vicinity of an interface cutting a cell at Xu,- 

In order to alleviate notation complexities we consider specifically cell 1 to be cut by 
an interface (wall) with notation illustrated in Figure [U at position x^, where a boundary 
condition is applied. We will consider the interpolation and solution on the left of the 
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interface, but the right is (under most conditions) simply the mirror image. We are seeking 
a solution of the potential on the whole grid: both sides of the interface. 



2.1.1 Robin Boundary Condition 

The boundary condition to be applied at we will take initially as being given by 

+ B(Pl + C = 0. 



(1) 



This form, sometimes called the "Robin" boundary condition, accommodates setting the 
value of 0, or 0', or some logarithmic gradient (j)'/4>. 

In bulk mesh branches (not affected by interfaces), the value of 0" applies only up to 
Xj+i/2. If an interface intersects the mesh between Xi and Xi+1/2, as illustrated in Fig. [H 
then 0" applies only up to that interface position, Xw, and 0-' must be evaluated using 
information from the boundary condition, not Xj+i, because Xj+i is not part of the solution 
to the left of the interface. If an interface intersects the mesh between Xj+1/2 and Xj+i, then 
again calculation of 0f cannot use Xj+i, and also 0" beyond this half-cell position cannot use 
standard differences and the values at Xi+1^2 because they are not part of the solution region; 
so we must make some other assumption. The most natural is to assume that the value 



(J)'- applies right up to the interface, 
derivative expression. 




Xu,, no matter its distance. Consequently, the second 
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must be considered to apply right up to x^. And in particular 

dxo + 2dxpi . , dxpi 
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We then determine 0^ and 0^ from a simultaneous solution of ([T]) and ([3]), and substitute 
into ([2]) to determine the second derivative, which after some algebra can be written: 
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{Adxpi + 2B) 



(5) 



where 

'^''p = ^'-'p^-^Adxpi + B) 

can be thought of as the equivalent mesh point distance corresponding to the boundary 
form. For B = (0^ specified) dXp is dxpi, while for A = (0^ specified) it is 2dxpi. 
It can easily be verified that the corresponding expressions for 0'/ are in agreement with 
intuition. This difference scheme is entirely equivalent to that of Shortley and Weller(2|, 
and subsequently Johansen and Colellaj^ (and others), when used with Dirichlet conditions 
{B = 0). These schemes are known to give V^0 to first order (truncation is 0{dx)) for 
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unequal differences!^, |5|], and hence automatically at boundaries. But because the errors 
are localized at the unequal difference points, the overall solution is generally second-order 
accurate jsj. 

Integrating, we obtain the equivalent interpolation: 

(f)'(x) = (f)[/2 + 01 (a; - X1/2), (6) 

and 

(j){x) = 01 + (p[/2ix - Xi) + (t)'[[{x - - {xi - Xil2f]/2, (7) 

which can be used to approximate the potential right up to the wall. 

Notice that this scheme uses a Taylor-expansion type interpolation, and appears to be the 
most accurate such scheme (for Laplace's equation or slowly varying second-derivative) that 
restricts itself to a three-point stencil. The three points are used to deduce three coefficients 
of the interpolation. Potential and its derivative are continuous except possibly at the wall. 



2.1.2 Media Boundary and Surface Charge Condition 

A different condition is to take the interface to be the boundary between two dielectric 
media with dielectric constant ei and 62. It is simple to generalize this condition to include 
a specified surface charge density cqct such that at the interface, the values in the left and 
right regions (1 and 2) are related by a jump condition 
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[Thus (J is the external charge in excess of any that arises from dielectric polarization.] In 
order to express this condition symmetrically, we need to extrapolate to get (p'^i using ([3]) 
but we also need to extrapolate from the right to get the gradient immediately on the other 
side of the interface 0^2- We do not wish to use an expression corresponding to ([3]) to the 
right of the boundary, because that would involve 03, and lead to a stencil that extended 
more than one leg from the node (1) under consideration. Therefore, we make the ansatz 
that the charge density (i.e. e0") has the same value in region 2 as in region 1. This is 
consistent with our presumption already made that 0" is uniform on one side in the vicinity 
of the boundary. The resulting boundary condition is 



ei0l,i - e20^„2 



which is 



dx, 
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ei0;'^ 
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Substituting this (0^, — 0i) into eq. ([2]) we find 

^„ _ (TdXp2/e2 + (02 - 0l) 

^ |_ dxpi + dxp2ei/e2 



£2- 



^2-0«, , ,f,dXp2 



dx. 



+ ei<; 



pi 



n)-T^ ei0i— + o- . 



dx 
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dXn 
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(9) 
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where 

dxp = dxpi + dxi/[l + {dxpie2)/ idxp2ei)]. (12) 

This gives the difference form and the gradient extrapolation through eqs. ([HI E])- It has the 
merit that when ei = €2 the differences reduce to those in the absence of an interface, and 
regardless of the ei, €2 values that the sum of the two dXp values adjacent to either side of 
the interface is always equal to 2xi. Naturally a gives an extra term representing a charge 
density. 



2.2 Difference Stencils 

In order to specify the difference scheme implementation for Poisson's equation 

= P, (13) 

it is helpful to have explicit formulas for the difference stencil. In the present approach the 
stencil involves only adjacent points. In A^^^ dimensions, this choice gives rise to a (2A^rf + 1)- 
point stencil. With reference to Fig. [H the difference equation at node (1) is written: 

(5: P5)0i - E ^505) = pi + E (14) 

s s s 

Here, S takes on the values corresponding to the adjacent nodes (negative, point 0, and 
positive, point 2, in the one-dimensional case being illustrated); Ps, Qs, and ts are coefficients 
that will be specified in a moment. Note that only the sums J^sPs and J2s^s need to be 
stored in addition to all Qs. Thus the storage per node required to specify the difference 
equation is 2Nd + 2 quantities. 

Along the dimension under consideration, the value of the diagonal coefficient is 

(15) 



{dXp + dXm)deff 

The (ixp^m-quantities can be considered to represent the distances from the node undercon- 
sideration to the control position in the positive dXp and negative dxm directions, and there is 
an extra effective distance d^ff. When there is no boundary, for example when the negative 
direction is between xi and xq, then the control distance is to the adjacent node: 

dxm = dxQ. (16) 

In the direction of the boundary, the other quantities are given in tabled! now to be explained. 

The "None" condition means a situation where there is no boundary. This row gives 
the standard coefficients of the difference scheme away from boundaries. The table specifies 
only the coefficients for the positive 6 (point 2 relative to point 1). The coefficients for the 
negative S (point relative to point 1) follow from refiectional symmetry. For example for the 
"None" case the negative side coefficients are given by eq. (fT5|) with dxm = dxi, dxp = dxo, 

deff = dXQ. 
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Boundary Condition 
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Table 1: Coefficients of the difference scheme. 



The Robin case, Aip + Bcf)' + C = 0, prescribes the condition that apphes on the left 
side of the interface position (i.e. it is what point 1 sees in the positive direction towards 
point 2). In this case, the value of potential at X2 is unused, and instead a constant r 
contribution is added. It acts like an additional charge density. Together the coefficients of 
this row are equivalent to equation (jl]). When there is an interface on the opposite side of 
the stencil center-point from the coefficients being evaluated, the value of dxm to be used 
is equal to dxp for that opposite direction. In other words, the coefficients for negative 5, 
in the stencil centered at xi, are obtained by taking dxp = defj = dxo, and Q = P, but 

"•^^ Adxpi+B "-^Pl- 

When a gradient boundary condition is applied using the Robin form with B ^ 0, then 
the boundary condition applied on the other side of the interface is usually just continuity 
of potential. If any other condition is applied, including the identical Robin condition, 
then a discontinuity in the potential will arise unless the exterior boundary conditions are 
consistently chosen. A discontinuity in potential gradient is usually physically acceptable 
(as a surface charge) but a potential discontinuity is usually not. Therefore a "Continuity" 
boundary condition is provided explicitly in the table to obtain the coefficients for point (2) 
relative to point (1) for the case where the Robin condition applies on the other side of the 
interface. That is, when point (2) sees the Robin condition in the direction of point (1) with 
the prescribed A, B, and C values. Naturally, when B = 0, this continuity condition reverts 
to the fixed-0 boundary condition: it becomes identical to the row above it. 

The "Media" row documents the stencils needed to implement the treatment explained 
in section [2.1.2[ 

3 Multiple Dimensions 

The one-dimensional problem just discussed is relatively straightforward. Extension to mul- 
tiple dimensions requires that a number of important complicating factors be considered. 

3.1 Difference Stencil and Boundary Conditions 

The difference scheme considered here for the Laplacian in dimensions is based upon the 
simplest stencil employing only the 2Nd adjacent points, giving a {2Nd+l)-pomt stencil. If 
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we consider a second-order Taylor expansion about a point that for notational convenience 
we take to be the origin of coordinates, it may be written: 

Nd Na Na 

= 01 + + + mi raf3XaX0, (17) 

a a a jjya 

where a refers to the dimension number and p, g, and r are the expansion's coefficients 
numbering A^^^, A^^^, and Nd{Nd — l)/2 respectively. The 2Nd differences of the (2Ai'rf+l)-point 
stencil permit the determination of only the pa and in addition to 0i. The r coefficients, 
which exist for A^^ > 1 in increasing numbers, cannot be determined. Fortunately, those 
cross terms do not contribute to the Laplacian, so that the multidimensional stencil can still 
determine the Laplacian to the same truncation order, using the 1-D difference scheme for 
each dimension. 

If the boundary conditions are of the Dirichlet type, specifying 0, then the generalization 
essentially straightforward. We use the B = version of eq (jlj) for second derivatives, and 
stencil contributions, in all coordinate directions. Their sum comprises the Laplacian. Such 
an approach in two-dimensions is frequently called the "Shortley-Weller" approximation, and 
is known to be second-order accurate, for essentially the same reasons as for 1-D. 

Gradient Boundary Conditions. If, however, we have a full boundary condition at the 
wall of the form 

A(f) + Bh.V(j) + C = 0, (18) 

where n is the normal unit-vector, then we need to translate this requirement into an equiv- 
alent for each of the dimensions, where wall intersections occur between nodes. It is the 
presence of the normal gradient d(j)/dn that is problematic, so Neumann conditions {A = 0) 
experience the same difficulty. Suppose the unit vector for (cartesian) dimension a is e^, 
and consider an elementary volume (half of a hypercuboid: in 3-D a tetrahedron, see Fig. 
|2]) defined by the planes normal to the dimension axes cut by the wall plane with normal 
n. If the area of the face of this volume defined by the wall plane is Q, then the area of 
the face perpendicular to axis-a is Qea.h. Consequently flux conservation demands that 
J2a^a-ii-'Va4> = n.V0. Therefore the flux-conservation implicit in the Laplacian is main- 
tained together with eq ( ITSll if in the a-dimension boundary condition, 

+ 5„V„0 + C = 0, (19) 

we take Ba = (-B/e^.n), (since Z]a(^a-n)^ = 1). This choice amounts to taking V0 = n|V0|, 
and ignores contributions to Va0 arising from any tangential component of the gradient. 
Such contributions might change the values of Va4> for individual directions a, but they sum 
to zero in the normal component. In so far as the control volume adjacent to the boundary 
can indeed be considered to be half a hypercuboid, the tangential gradient does not contribute 
to the divergence of the flux, because at the interface its normal component is zero and at 
the lattice cell-boundaries the flux is obtained directly from the flnite differences there, 
regardless of what is the flux tangential to the oblique interface. Thus, for multidimensional 
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Figure 2: Interface intersecting the lattice legs forming a tetrahedron in 3 dimensions. 



evaluation of the Laplacian, in the vicinity of an interface, one can use for each dimension 
the one-dimensional stencil derived in section |2] but with the boundary condition coefficient 
S„ = (5/e,.n). 

By the same argument, when a surface charge density a is present, which gives rise to a 
gradient discontinuity, the fraction of that charge density to be attributed to dimension a 
is e^.n. Then, associated with an interface area ^2, dimension a possesses a projected area 
f^eQ.fi. When each dimension area is multiplied by its charge density ea.ficr and the resultant 
summed over all dimensions, the total surface charge Via results. A lattice leg along the a 
dimension intersecting the surface must use the 1-dimensional difference scheme of eq. (fTTj) . 
but with charge density eQ,.no". 

These choices are the appropriate ones for a difference scheme with a (2A^rf-|-l)-point 
stencil in iV^- dimensions. With the standard expression for the second derivative along each 
dimension or its generalizations, eqs. (jlj) or (II ip in hand, (as summarized in Table 1) one 
can then solve the resulting matrix equation by standard methods such as SOR or conjugate 
gradient minimization, to find the potential on the mesh. 

3.2 Multidimensional Interpolation 

More complex mathematical considerations and by far the more difficult coding challenges 
lie in the back-interpolation of the solution to arbitrary position, rather than the solution of 
the difference equations. 

Gradient interpolation. For the dynamics of particles moving in the field, the gradient 
of the potential is the most important. Gradient interpolation is done most naturally using 
the information embedded in the difference stencils. For gradients in direction a at positions 
along a-lattice-legs, i.e. at coordinate values in the other dimensions corresponding to a 
node, the 1-D interpolation is entirely natural and straightforward. Eq. ([6]) with eq. (jlj) at 
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bulk nodes and Dirichlet boundaries, or with eq. ( ITT]) at Media and surface- charge interfaces 
can be used. 

However, at interfaces where the Robin (with i? 7^ 0) or any normal-gradient condition is 
applied, interpolation that sets the tangential gradient to zero is insufficient. It is sufficient for 
evaluating the Laplacian to ignore tangential gradient but not for the field interpolation. We 
must estimate the tangential gradients (A^^^ — 1 quantities). We wish to do this for convenience 
and compactness without extending our consideration beyond the standard {2Na + 1) stencil 
node values. We therefore have to drop cross-gradients from the Taylor expansion as we 
do in the bulk scheme. That means we can ignore the normal derivative of the tangential 
gradients, and determine the gradients on the opposite side of the central node from the 
interface intersections. That is (see Fig. |2]) from the gradients at positions Xa,i/2, taking the 
tangential field to be simply constant for the cell under consideration. The estimate of the 
tangential gradient is therefore 

Vt(f) = V(f)\i/2 - (n.V(/)|i/2).n (20) 

(where an appropriate interpolation for fi is used if curvature is present). Then when the 
interface boundary condition is the specification of the gradient at each of the face points 
(dxpi) is 

V(^U = + = V0|i/2 - (n.V0|i/2) .n + . (21) 



_ (22) 

Substituting for from the Robin condition we get 

- - -n.{A/B)^^ - ( C/B + E^,^^^ - 1 ■ (23) 




aw 



This is now in the form of a Robin condition in the single (a) dimension, but with new Robin 
coefficients 

A'^A, B' = B/n„, C' = C + B(^n,tl^.^zM . (24) 

These primed coefficients can be substituted into our prior expression for the second deriva- 
tive, eq dl]), to obtain: 

C + A<l>, + A(0i - <l>^o)'^ + B T.f n,^l 



AdXapl + B/Ua dXaO + dXp ' ^ ^ 

where dxp = dxapi{Anadxapi + 25)/ {Auadxapi + B). It is helpful to note that in the limits 
-B — !■ 0, or — > 1 this expression is identical to ([1]). For other parameters, the most 
important extra distinction is the inclusion of the B term in the numerator, which involves 
values (ban off the a-lattice-line. 
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Figure 3: Gradient-interpolation geometry within an elementary lattice box. 

Transverse interpolation Referring to Fig [3] which shows an example in 3 dimensions, 
we use as base the node closest to the point, and the "box" consisting of the elementary 
cuboid whose vertices are adjacent mesh nodes, in which the point lies. 

For the gradient along the a-direction, we use the quadratic interpolation, eq ([6]), along 
each of the 2^"*"^ lattice legs (i.e. box edges) aligned along the a-direction, to obtain the 
gradients at interpolated points where the a-coordinate fraction of the lattice 

spacing from the base, and the other coordinates correspond to nodal values. Then we use 
multilinear interpolation between those values in the orthogonal (A''rf-l)-dimension subspace 
to approximate the gradient at the point. Multilinear interpolation can be considered to be 
defined inductively via 

(j)'{Xa, XfS, X^, --XNa) = fv^f^'i^a, Xp, X^o, ■■■Xm) + (1 " fv)(P'i^a, XfS, Xr,l, ...Xm) (26) 

for all rj ^ a. It has the important merit of giving continuity at cell boundaries but requiring 
only localized box information. This scheme provides gradient (field) estimates that are 
continuous and have at least first order accuracy. 

If at the first, quadratic, stage we find a node from which we would do interpolation, e.g. 
(xqo, x^o), to be outside the point's solution region, i.e. on the other side of an interface, 
then we examine the other end of the lattice leg under consideration (xqi, X/31, x-yo)- If that 
end lies in the point's region, then we extrapolate from that instead. If that second end does 
not lie in the point's region either, then the lattice leg is entirely outside the interpolation 
region. In that case, an "over-extrapolation" is sought by looking one more node in either 
direction. If a node is found that is in the region, then (over-)extrapolation from that node 
to the lattice position under consideration is adopted. Such a value is considerably more 
uncertain than values corresponding to points that are in the region or directly extrapolated 
from the region. Over-extrapolated values are given a weight w equal to one minus the 
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fractional position relative to the lattice end adjacent to the region. That is, if position 
Xao is in the region, and Xai, Xa2 not, then the weight at x such that Xai < x < Xa2 is 
w{x) = 1 — (x — Xai)/{xa2 — x^i)- If there is no extrapolation available, then the weight is 
set to zero. 

When one or more of the weights in a multilinear interpolation is not unity, a weighted 
multilinear interpolation is used. It can be expressed as a sum over all the 2'^^''"^^ dimensional 
corners of the box within which the interpolation is being done, denoted by the vector index 
whose components have the values = or 1, where 1 < rj < N^ — l denotes the dimension. 



Thus 

0'= ^ Mh) where A{ir,) = w{ir,)Y[il - tr, + U2tr, - l))- (27) 



When all w = 1, this expression is equivalent to eq ( 126|) . 

For concave interfaces, it is theoretically possible for a particle to be in a box which has 
no vertices in its region. If that happens, a fall-back is needed, the natural thing is to choose 
as base node the nearest node that is in its region. This problem can usually be avoided by 
choice of mesh in relation to objects. 

Potential interpolation. Because we do not directly differentiate the potential, using 
instead the gradient interpolation just explained, satisfactory values for the potential are 
obtained by multilinear interpolation, and do not require higher order interpolation schemes. 
Multilinear interpolation is unproblematic in the "bulk" where all the box vertices are valid 
points in the relevant region of the grid, and no interface crossings occur on its edges. The 
problematic issue is how to extrapolate to the interface, when box vertices are absent from 
the valid region because they are the other side of an interface. After some experimentation, 
the following approach was found to work well using modest computational resources. First, 
all the vertices of the box are examined. If any do not lie in the region, they must be filled in 
with extrapolated values. The extrapolated values are obtained by calculating (during this 
examination) the centroid, the mean value, and the mean gradients of the "valid" vertices 
of the box (ones that lie in the region). Additional steps must be taken for any dimension 
in which the gradient is undetermined, which happens when all the valid vertices have the 
same coordinate in that dimension. Such a situation occurs when an entire face of the 
box is outside the valid region. The occurence is not all that unusual. In such a situation 
of undetermined gradient, we instead determine the gradient in that dimension, using the 
gradient interpolation code, evaluated at the point. The values that are filled in to the 
invalid vertices are then the mean cell value plus the dot product of the gradient and the 
vector to the vertex from the centroid. Multilinear interpolation can then proceed as usual 
on the entire box including the filled-in vertices. 
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Figure 4: Representation of spherical object on 16^ (a) and 32^ (b) grids. 

4 Example Accuracy Tests 

Here we give examples of the typical accuracy of the difference solution and interpolation 
schemes. The case illustrated is of a sphere of fixed potential 10, radius 1, inside a second 
sphere of radius 5 at which the potential logarithmic radial gradient is -1. It is solved on a 
mesh over the region —5 < x,y, z < 5. The inner sphere wireframe is shown in Fig. Hlfor 
two grids, 32^ and 16'^. Also shown in this figure are the lattice legs which are intercepted 
by the sphere. The wireframe nodes are those intercepts. 

The solution of this problem: = 0, is of course (j) = 0(l)/r. In Fig. [5] the values at 




-6.0 -4.0 -2.0 0.0 2.0 4.0 6.0 
X mesh (at half mesh in y, z) 



Figure 5: Solution of the difference equation, 
the nodes in the plane of constant y, z (through the center) are compared with this analytic 
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solution. Inside the inner sphere the potential is constant, dictated by continuity at the 
interface, as would occur for a conducting sphere. (Outside the outer sphere the potential is 
set to zero.) 

The values at the mesh nodes tell only the accuracy of the solution representation at 
the nodes. Of more interest is the overall accuracy of the solution when interpolated to 
positions away from the nodes. Fig. [6] illustrates this by performing the interpolations 
previously described along a line of constant direction in the z = plane whose angle is 
0.1 radians (5.7°) to the x-axis. The lines labelled "getpotential" and "getfield" are the 




0.0 1.0 2.0 3.0 4.0 5.0 0.0 1.0 2.0 3.0 4.0 5.0 



Figure 6: Comparison of different interpolation schemes for 32^ mesh. 

quadratic interpolations. For comparison, the analytic solution is plotted and two other 
types of approximation. That called "nearest" simply takes the potential to be the value 
at the nearest node of the mesh. The approximation called "simple" is arrived at by a 
simple linear interpolation of (p between nodes, and for "simplefield" a linear interpolation 
of V0 based on taking differences of the (p nodal values. The "simple" and V(/> estimates 
therefore ignore the extra information embedded in the boundary conditions. We notice 
that although the (p interpolation in the bulk region is very good for all approaches, near 
the boundaries there are big errors for the simple and nearest approximations. For the 
getpotential interpolation, the errors are an order of magnitude smaller, having a maximum 
error of 0.145 along this line, occurring at a position just outside the inner sphere. This 
value is typical. (The "region/10 line" simply indicates the different regions of the solution.) 

The field (V0) interpolation shows even larger errors in the simple estimate, because 
linearly interpolating it from differences in the vicinity of its discontinuity at r = 1 produces 
a piecewise linear continuous behavior that is far from the analytic form. The quadratic 
interpolation using the boundary information, by contrast, accommodates the discontinuity 
in field and obtains very respectable approximations even in the vicinity of the boundary. 
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The maximum field error varies substantially with the angle of the path. In this case it is 
about 0.8 in radial field component and less in tangential field. 

Figure [7] documents the spatial dependence of the error in the x-direction field in a plane 
y =constant ~ 0). The contours of error are spaced by 0.2. The positions of the nodes 




-1.0 -0.5 0.0 0.5 1.0 
axis-3 



Figure 7: Error in field component in direction-1 within a plane of constant y (axis-2). 
Contour spacing 0.2. (The maximum of |V0| is 10; so contours are 2% relative error.) 

are marked by squares. What we see is that the field error is less than 0.2 over most of the 
domain but locally near the field discontinuity at r = 1 larger errors occur that depend upon 
position round the boundary. 

The scale length of the field and potential variation immediately outside r = 1 is 1; the 
mesh spacing is about 0.3. It is therefore evident that we are obtaining highly favorable 
accuracy, typically 10% maximum local error and 1% standard deviation, in field even with 
rather coarse mesh resolution (0.3) of the scale-length. Moreover using the boundary infor- 
mation we completely avoid pollution of the field in one region by values in another where 
it is different by a discontinuity. 



Nodes 


16^ 


32^ 


483 


643 


1283 


Spacing dx/vs 


0.63 


0.31 


0.21 


0.156 


0.078 




0.020 


0.0031 


0.0014 


0.0008 


0.0002 


^-Ejmax / EjYidx 


0.167 


0.091 


0.068 


0.050 


0.022 


SD{AE)/Emax 


0.020 


0.0089 


0.0060 


0.0022 


0.0007 



Table 2: Nodal-potential and field-interpolation relative error scaling with mesh size. 

Table [2] show the maximum nodal (i.e. uninterpolated) potential error A(f)max/4>max, and 
the maximum local (interpolated) field error AEmax/Emax and standard deviation of the 
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error SD{AE)/E^ax evaluated over a central cube of side length 3. These errors decrease 
with the mesh spacing, as the number of nodes is increased. The nodal potential error 
decreases quadratically with mesh spacing as expected for this second order scheme. The 
maximum interpolated field error decreases only linearly, but the error becomes more and 
more localized so that the standard deviation decreases approximately quadratically. By 
comparison, the unsatisfactory "simple" interpolation (not shown in the table) gives local 
field errors whose maximum magnitude is independent of mesh spacing. The "getfield" 
interpolation asymptotic scalings arc the same as a "nearest" value interpolation (using 
the nearest value of the gradient that is based on values only from nodes on the same 
side of any boundary as the point being interpolated). This is presumably because there 
are situations near boundaries where the multilinear interpolation provides no tangential- 
field derivative-information, omitting nodes that lie across a boundary. The magnitude of 
getfield's error, however, is smaller by a substantial factor (perhaps 3 or 4) than a nearest 
value approximation, and the volume of its region of bad approximation is smaller. 

In summary, using a compact difference stencil, Poisson solution and interpolation of 
the electric field is implemented with excellent accuracy and convergence, suitable for use in 
particle in cell codes with purely cartesian mesh but essentially arbitrary oblique boundaries. 
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